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A nonequilibrium distribution function of microscopic thermal current is studied by a direct 
numerical simulation in a thermal conducting steady state of particle systems. Two character- 
istic temperatures of the thermal current are investigated on the basis of the distribution. It is 
confirmed that the temperature depends on the current direction; Parallel temperature to the 
heat-flux is higher than antiparallel one. The difference between the parallel temperature and 
the antiparallel one is proportional to a macroscopic temperature gradient. 
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The characterization of nonequilibrium steady state 
(NESS) is one of the important issues of nonequilib- 
rium statistical mechanics. For a linear nonequilibrium 
regime, several methods describe the nature of NESS. 
The linear response theory gives a method of calculat- 
ing response coefficients from equilibrium fluctuations. 
An approach using the Boltzmann equation enables us 
to obtain response coefficients through a nonequilibrium 
distribution. However, in general situations such as that 
in a nonlinear nonequilibrium regime, we have no clear 
nor general methods for such a description. 

In the following, NESS is analyzed on using a distri- 
bution function, because a distribution function implies 
direct information on NESS. In particular, we investi- 
gate a microscopic distribution of the energy current of 
particle systems in a steady thermal-transport state. A 
thermal-transport system is a simple and well-studied 
example of nonequilibrium steady states. 1 Microscopic 
energy current is carried by a single particle in a particle 
system. If no net flow of particles exists in the system, 
the energy current corresponds to microscopic thermal 
current. Gathering microscopic thermal current in an ap- 
propriate cross section, we obtain macroscopic thermal 
current. From the viewpoint of nonequilibrium statisti- 
cal mechanics, macroscopic thermal current has less in- 
formation than microscopic thermal current. Generally, 
its distribution is Gaussian because of the central limit 
theorem and has only two parameters: average and vari- 
ance. 

Let us consider an equilibrium case of the microscopic 
distribution before a nonequilibrium case. For a particle 
system in three-dimensional space, an ^-component of 
the kinetic part of the thermal current is defined as 2 
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with a mass m and a momentum p = (p x ,p y ,p z ) ■ The 
distribution of is obtained by taking the thermal av- 
erage of S(j — ijj) with a Maxwellian distribution of tem- 
perature T. (Hereafter, we take the Boltzmann constant 
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to be unity.) It is given by 



P eq (j) = (S(j 



iS)) eq = {i^ (* 



with a normalization <i?'P eq (j) = 1, where fi is a 
temperature-scaled mass parameter defined by fi = m/ 
T 3 , which has the dimension of the inverse square of ther- 
mal current, and a dimensionless thermal current with 
a fractional power z = (fi/2) 1 ^ 3 \j\ 2 ^ 3 . Ei(z) is the En- 
function with n = I defined by 3 



Ei (2) 
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which is related to the exponential integral function 
Ei(z) = — Ei(— z). In this letter, we treat only the ki- 
netic part, although a distribution of the potential part 
of the thermal current can also be evaluated. 4 The equi- 
librium distribution eq. (2) has a log divergence, — lnz, 
near j ~ and an exponential tail with an algebraic 
correction, ~ e~ z / z (z — > oo). 

Let us consider a nonequilibrium distribution of 
For a small thermal current j ~ 0, it should behave like 
an equilibrium distribution with a local equilibrium tem- 
perature, that is, it has a log divergence, because the 
small current is contributed by low-energy particles that 
are well thermalized locally. High-energy particles con- 
tribute to the negative and positive tails of the distri- 
bution. They come from far places from the observation 
point without scattering by other particles. Therefore, 
these high-energy particles have information on these 
places. When these places are separated by a scale larger 
than the mean free path, the tails of the distribution of 
the thermal current may have some information on local 
equilibria. If these local equilibria are characterized by 
two typical equilibria that are separated by the scale of 
the mean free path, it means that the apparent tempera- 
tures of the tails are different from the local temperature 
of the observation point, and that the two temperatures 
are related to a thermal gradient and a scale of the mean 
free path. If the temperature gradient dT/dx is negative, 
a large positive x component of thermal current has in- 
formation on the high-temperature local equilibrium and 
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Fig. 1. Geometry of nonequilibrium particle simulation: The 
component of the kinetic part of thermal current is measured 
in the cross section of the system with the width a located at 
the center of the box. 

a large negative x component has information on the low- 
temperature local equilibrium. 

The above picture predicts that the tails of the dis- 
tribution of the current can be described by two other 
local equilibrium distributions with different tempera- 
tures. For an ideal gas system, this is trivial. In this 
letter, we study this picture for dilute Hertzian-sphere 
and Lennard- Jones particle systems by direct numerical 
simulation. 

In order to confirm the picture of the tail structure, we 
perform a nonequilibrium particle dynamics simulation. 
The geometry of the system is shown in Fig. 1. Particles 
are confined in a rectangular parallelepiped box, whose 
size is denoted the L x x L y x L z . Heat bath regions are 
attached to both sides of the x-direction with a width Lb, 
where the kinetic temperature of particles is controlled 
by a Nose-Hoover thermostat 5-7 with different temper- 
atures (Th,Ti). The boundary conditions are as follows: 
In the x-direction, we impose appropriate elastic wall po- 
tentials on both ends; in the y- and z-directions, periodic 
boundary conditions are imposed. The interaction poten- 
tial for the inter-center distance r of particles is taken to 
be a hard Hertzian potential, 8 <j>(r) = Y\a— r| 5 / 2 8((7— r), 
where Y is proportional to Young's modulus, which is 
taken to be Y = 10000e/cr 5 / 2 with an energy unit e, a is 
the characteristic length of the particles corresponding to 
a diameter of Hertzian particles, and O is the Heaviside 
step function. Empirically, this potential can reproduce 
the properties of a hard-sphere system. The dynamics 
of the particle is governed by Newtonian dynamics in a 
bulk system. Using heat baths with different tempera- 
tures, we realize a nonequilibrium steady thermal con- 
ducting state in the simulation. From previous studies, 
this setup enables a normal thermal conduction, which 
is characterized by the Fourier law (linear profile of local 
temperature) and a finite transport coefficient. 9-12 

In the following simulations, we take the following pa- 
rameters: L x = L y = 16cr, L z = 120cr, and Lb = 16a. 
The number of particles N = 3072. We take three 
different pairs of heat bath temperatures: (Th,T{) = 
(lOe, e), (20e, e), and (40e, e). Discarding the initial re- 
laxation steps, steady states are realized. Their local 
number density n at the center is about n ~ 0.1 /cr 3 in 
all cases, and their local kinematic temperature T at the 
center are T ~ 6e, lie, and 23e, respectively. The thermal 
gradients of actual steady states near the center are also 
calculated as dT/dx ~ 0.06e/cr, 0.12e/<r, 0.23e/er, respec- 
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Fig. 2. Distribution of j^, which is a parallel component of the 
thermal gradient: The horizontal axis represents sgn(j)z and the 
vertical axis represents ln(P(j)/ y/Jl), where z and fi are identical 
to definitions used in the equilibrium distribution cq. (2). Three 
different thermal gradient cases are superposed. 




Fig. 3. Negative tails are flipped at z = and superposed over 
positive tails: Gray scattered dots are original data and black 
lines are averaged over several bins for each data point of three 
different gradients. Clearly, we observe two branches for negative 
and positive currents. 



tively. Note that the log-thermal gradients of the actual 
states are dlnT/dx ~ O.Oler -1 in all cases. 

The distribution of j£ is observed in the cross sec- 
tion of the system with a width a located at the center 
[60cr, 61er). The results for about 3 x 10 6 samplings per 
case are shown in Fig. 2. We recognize that the distri- 
butions of jj| with three different thermal gradients are 
well scaled with sgn(j)z and P{j)/^JJi from the figure. 
This scaling collapse of the distribution for three different 
thermal gradients is consistent with the analysis based on 
the Boltzmann equation. In the Boltzmann analysis, the 
thermal gradient always appears in the form, dhiT/dx, 
instead of dT/dx. 13 In the present simulation, the log- 
thermal gradients are 0.01 in all cases. Thus the scaling 
collapse means that the simulation data have the same 
dependency on the thermal gradient as the Boltzmann 
equation. The simulation also confirms that the distri- 
butions of the other two components are identical to the 
equilibrium distribution P eq (j). 
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Fig. 4. Scaling plot of the negative and positive tails: Gray scat- 
tered dots are scaled original data and black lines are scaled 
branches for three different thermal gradients. Six curves are 
well collapsed into a single curve. 



The tails of the distribution in Fig. 2 are slightly 
skewed. To investigate the form of the tail in detail, we 
replot the data by flipping the negative tail at z = and 
superposing the positive tail with statistical smoothing 
over several bins. The result is shown in Fig. 3. Plots 
from three different thermal gradients clearly show the 
same curve again. In addition, branches of positive and 
negative tails are clearly distinguished. 

If the conjecture that the tails of the distribution are 
identical to equilibrium distributions with different tem- 
peratures is true, we can expect another type of scaling 
from the form of the equilibrium distribution eq. (2): 



p(j) 3/2 n?) 



a± 
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where a± is a rescaling factor of local temperature. For 
the positive tail, its temperature is given by T> = a + T 
and, for the negative tail, its temperature is given by 
T< = a_T. These temperatures are called "scaling" tem- 
peratures. The best fit of the above scaling is given by 
a± = 1 ± 0.04 and the result is shown in Fig. 4. This 
scaling factor is also related to the thermal gradient dT/ 
dx and the characteristic length I as 



a±T 



dT, 

T± 7T l 
dx 



(5) 



a± = 1 ± 0.04 and dlnT/dx = O.Oler" 1 give the char- 
acteristic length I as I = 4<r. The mean free path If of 
the present situation is calculated as If = l/y^nna 2 ~ 
2.25cr. Thus, the picture of nonequilibrium tails are con- 
sistent with the computational result. 

Here, we use another definition of tail temperature. 
Calculation of the new defined temperature is easier than 
calculation of the scaling temperature: Local-equilibrium 
kinetic temperature T is calculated by 

1 ' (6) 
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where (. . . ) ne ss means that the average is calculated in 
the local section in NESS. Decomposing this average into 
two parts with positive and negative p x ensembles, we 
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Fig. 5. Difference in effective kinetic temperature between both 
tails normalized by the local kinetic temperature plotted against 
the local temperature gradient for the Lennard-Jones particle 
system. 



can define the new temperature T e g (T g) for the positive 
(negative) tail: 

^ = ^(P a >a>°«odis = ^ 

From the definition, the following relation exists: T = 
(T e g + T c ^)/2. This temperature is a sort of "effective" 
temperature because of its difference from the scaling 
temperature. The difference comes from contributions 
near j ~ 0. For the present results, T^/T and T^/T, 
which correspond to the alpha parameters of the scaling 
temperature, are 1 ± 0.0235(5), 1 ± 0.0284(4), and 1 ± 
0.0335(5), for dT/dx ~ 0.06e/cr,~ Q.12e/a, and 0.23e/ 
a, respectively. With increasing thermal gradient, the 
difference between the two temperatures becomes much 
smaller. 

For the Lennard-Jones particle system, a similar result 
is also observed. The Lennard-Jones interaction is taken 
to be as follows: 



(8) 



with a potential cutoff length 3r c . Simulations with the 
same setups as the Hertzian particle system are con- 
ducted with the following parameters: L x = L y = L z = 
15r c , and — 2r c . The number of particles is N — 
843. Distribution is measured at the center [7r c ,8r c ). 
The average number density is n ~ 0.249/r^. Several 
temperature gradients, (T^,T;) = (3.5e, 3.3e) (3.6e, 3.2e) 
(3.8e, 3.0e) (4.0e, 2.8e) (4.2e, 2.6e) (5.0e, 1.8e), are exam- 
ined. These parameters are chosen as the local kinetic 
temperature in the observation section becomes about 
3.4e. From these parameters, the Lennard-Jones system 
is in a supercritical fluid phase. The simulation gives us 
a skewed nonequilibrium distribution of The temper- 
atures of the tails also differ from the local equilibrium 
temperature. We can calculate the scaling and kinetic 
temperatures. Again, there are differences between both 
temperatures. For the Lennard-Jones system, we plot the 
difference between T c ^ and T c ^ normalized by T against 
several temperature gradients in Fig. 5. In the small- 
gradient region, a linear relation is observed. It seems 
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that linearity is broken in a rather large thermal gradi- 
ent. At this thermal gradient, the linearity of thermal 
transport is also broken. 

Such kind of tail property in which the negative and 
positive tails obey the equilibrium distributions with dif- 
ferent temperatures asymptotically are not special prop- 
erties of particle systems. Recently, for the quantum ther- 
mal transport system of the harmonic chain, it is shown 
that the distribution function of heat has exponential 
tails with different temperatures. 14 For the stochastic lat- 
tice thermal-conduction model proposed by Giardina et 
al. 15 an exact nonequilibrium distribution of microscopic 
thermal current through a bond connecting site i and site 
i + 1 can be calculated. The result is 
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where Tj and T i+ i are the local temperatures of subse- 
quent sites. The equilibrium distribution of the Giardina 
et al.'s model with a temperature T is expressed as 
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The nonequilibrium distribution has the same tail prop- 
erty, because of the asymptotic behavior of the modified 
Bessel function K (x) ~ e~ x j\fx, as x — > oo: 
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Therefore, the behavior of the tail, in which the nega- 
tive and positive tails obey the equilibrium distributions 
with different temperatures asymptotically, is a univer- 
sal property of the nonequilibrium steady thermal trans- 
porting system, at least, in the linear nonequilibrium 
regime, although the asymptotic form itself depends on 
the dimensionality of the momenta and details of the 
models. 

In this study, we have investigated the microscopic 
distribution of the kinetic part of heat current by di- 
rect nonequilibrium numerical simulations of particle 
systems. The result shows that the positive and negative 
tails of the distribution for the parallel component to the 
thermal gradient have an asymptotic form of the equi- 
librium distribution with different temperatures. These 



temperatures differ from the local equilibrium temper- 
ature. The temperature of the positive tail, which cor- 
responds to the normal current direction, is larger than 
that of the negative tail. In addition, we have found that 
these temperatures are identical to those of the places 
separated by a distance several times the mean free path 
scale. Two definitions of the temperature of the tail have 
been studied: One is the scaling temperature, which is 
accurate but difficult to obtain, and the other is the ef- 
fective kinematic temperature, which is less accurate but 
easy to calculate. Both temperatures can capture of the 
tail structure well. 

The tail structure is not a special property of the par- 
ticle system, but it may be a universal property of a nor- 
mal thermal transport. For a nonlinear regime of ther- 
mal transport, the same tail structure of has already 
been observed in the Lennard- Jones particle system. In- 
vestigations of the potential part of thermal current that 
consists of potential advection and work terms are now 
in progress. 
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